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We use fully nonlinear numerical relativity techniques to study high energy head-on collision of 
nonspinning, equal-mass black holes to estimate the maximum gravitational radiation emitted by 
these systems. Our simulations include improvements in the construction of initial data, subsequent 
full numerical evolutions, and the computation of waveforms at infinity. The new initial data 
significantly reduces the spurious radiation content, allowing for initial speeds much closer to the 
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obtained by thermodynamic arguments (0.134) and by previous numerical estimates (0.14 ± 0.03). 
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I. INTRODUCTION 

The study of the high energy collision of two black 
holes is of interest from both the theoretical point of view, 
to understand gravity in its most extreme regime, and ex¬ 
perimentally, since increasingly high energy particle col¬ 
lisions could eventually have a non-negligible probability 
for generating black hole pairs (see Ref. [1] for a review). 

The production of gravitational waves and the proper¬ 
ties of the final remnant after the collision of two black 
holes has been the subject of theoretical study for over 
half a century, with notable results such as the area theo¬ 
rems by Hawking and Penrose w and their application 
to bounds on the energy radiated via gravitational waves. 
For instance, they find an upper bound for the maximum 
energy radiated from a head-on collision of nonspinning 
black holes of 29% of the total mass. 

More detailed estimates of the radiated energy have 
been computed by applying perturbation theory [3] to 
the collision of ultrarelativistic black holes represented by 
shock waves [S]- Those computations reduce the above 
bound to 25% (when only including first-order correc¬ 
tions) and to 16.4% (when second-order corrections are 
included. A ZI—dimensional generalization of the first- 
order computation [S] found that the proportion of en¬ 
ergy radiated to the initial mass scales as 1/2 — 1/D. 

Fully nonlinear numerical simulations of such collisions 
are now possible thanks to the breakthroughs in numer¬ 
ical relativity OE]. The first full numerical study of the 
head-on collision of black holes uni found a maximum 
efficiency of 14 ± 3%. Those studies have been extended 
to grazing collisions leading to an estimate of 35% 
for the maximum energy radiated at a critical impact pa¬ 
rameter. Further studies including boson stars [12] , fluid 
stars mm], black hole spins m and unequal mass bi¬ 
naries [16] show that, at high energies, the structure (i.e. 
matter, spins and mass ratios) of the holes tends to be 
irrelevant for the collision outcomes. 


The latest analytical computations of the energy ra¬ 
diated by the head-on collision of two, equal mass, non¬ 
spinning black holes include an estimate of 13.4% based 
on black hole thermodynamics arguments [17] and 17% 
based on a multipolar analysis of the zero-frequency-limit 
(ZFL) approach [15] . 

In this paper we revisit the full numerical head-on 
computation incorporating new techniques that improve 
the accuracy of the simulations. These techniques in¬ 
clude new initial data with reduced spurious radiation 
content m, improved extraction techniques with sec¬ 
ond order perturbative extrapolation [50] , and the use of 
new gauges [21] and evolution systems [22] in the moving 
puncture approach [8]. 

We use the following standard conventions throughout 
this paper. In all cases we use geometric units where 
G = 1 and c — 1. Latin letters (i, j, ■■■) represent 
spatial indices. Spatial 3-metrics are denoted by and 
extrinsic curvatures by Kij. The trace-free part of the 
extrinsic curvature is denoted by Aij. A tilde indicates 
a conformally related quantity. Thus jij = ^^^d 

Aij = 'tp~‘^Aij, where ip is some conformal factor. We 
denote the covariant derivative associated with jij by Di 
and the covariant derivative associated with by Di. 
A lapse function is denoted by a, while a shift vector by 
(3\ 


II. NUMERICAL TECHNIQUES 
A. Initial Data 

We use an extended version TwoPunctures [ 23 ] 
thorn to generate puncture initial data m for boosted 
black hole binary simulations. In the conformal 
transverse-traceless (CTT) formalism [24ll27j . the con¬ 
straints on the initial spatial hypersurface Eq become a 
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set of elliptic differential equations for the conformal fac¬ 
tor and potential vector (see Eq. © below) through a 
conformal transformation 


The trace-free part of the extrinsic curvature is split 
into background terms and a longitudinal correction 
term obtained from a vector b^. Here 


lij = ■ 


Mij — 


i+) 


■ a: 


(-) 


( 8 ) 


We call ^ij the conformally related metric tensor. All 
objects with a tilde are associated with 7 ^-. 

As in Ref. m. to calculate the spatial metric and ex¬ 
trinsic curvature associated with a boosted black hole of 
mass m and arbitrary linear 3-momentum P*, we Lorentz 
boost the 4-dimensional Schwarzschild line element in 
isotropic Cartesian coordinates. We then extract from 
the transformed metric the spatial metric 7 *^-, the lapse 
function a*, and the shift vector /3* (a super/subscript * 
indicates that this is a single black hole quantity). We 
then obtain the extrinsic curvature K*^ on Eq using the 
evolution equation for the spatial metric 

Kh = ^{D:p*+Dp:-dtn:j) ■ 

CTT separates this into trace and trace-free parts 



where K* = K*j. For the conformal factor, we make 
the standard choice 




( 1 ) 


where r is the unboosted isotropic radius. 

For example, if the boosted coordinates are given by 


t' = "ft + jvy , 
x' = X , 

y' = iy + ivt , 


( 2 ) 

(3) 

(4) 

(5) 


then the conformal spatial line element on Eq (defined 
by t' = const) is given by 


de'^ =dx'‘^+j 


16{m — 2r)^r^v'^ 
(m -|- 2 r)® 


1 - 


dy'^+dz'^ , ( 6 ) 


where v is the magnitude of the local velocity vector 


P* 

P = , (7) 

-h pi Pj 

(here the boost is along the y-axis), 7 = (1 — 
and r = ^x'^ + y"^ + z'^ = x'"^ + — vt'Y + 

Our black hole binary initial data is constructed using 
a superposition of metric and extrinsic curvature terms 
derived from the above expressions. To distinguish con¬ 
tributions for the two black holes, we replace the * su¬ 
per/subscript above with a -I- or —. 


where and A,-■ ^ are the trace-free part of the con- 
formal extrinsic curvature of a single boosted black holes 
located at r = r+ and f = r-. Note that the trace- 
free part of the single boosted black hole extrinsic cur¬ 
vature will have a small trace with respect to a metric 
constructed by superimposing two different background 
metrics. We remove this extra trace term prior to solving 
the initial data equations, i.e., Mij —>■ 

(where 7 ^^ is the superimposed background metric). The 
complete trace-free part of the extrinsic curvature for the 
superimposed spacetime is given by 

Aij = Mij + —{i^b)ij , (9) 

where a = ip^d and {U)ij = Dibj + Djbi — ‘^^ijDkb^ 
is the longitudinal vector gradient. As part of the freely 
specifiable parameters we set a = 1 . 

In the puncture approach, we write the conformal fac¬ 
tor as singular parts plus a finite correction, u, 


■0 =' 0 (+) + V'(-) - 1 + w , ( 10 ) 

where ^-re the conformal factors Q associated with 
the individual, isolated black holes located at positions 
labeled as (-I-) and (—), with spatial metric tensors 
Given these choices, the Hamiltonian and momentum 
constraints become equations for the correction functions 
u and 6 * 
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+ D^ (V'(+)+V'(-)) = 0 , 

( 11 a) 


A DjM"^ - = 0 , 

' (lib) 


where A]l 6 * = Dj{hb)d is the vector Laplacian and R is 
the scalar curvature associated with 7 ^-. The solutions 
are required to obey Dirichlet conditions at infinity 

lim u = 0 and lim &* = 0 . 

r—yoo r—>-oo 

In order to deal with the puncture singularities, we in¬ 
troduce attenuation functions to both modify the back¬ 
ground metric and mean curvature, as well as modify 
the singular source terms inside the horizons themselves. 
The first type of attenuation, which is consistent with 
the constraints everywhere, is used in the superposition 
of the background conformal metrics and has the form, 

7ij = Sij + /(+) ( 7 !/^ - Sij'j + /(_) ( 7 !^^ - Sij'j , 
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where 


B. Convergence of Initial Data 


/(±) = 1 - , 


and r(j-) is the coordinate distance from a field point to 
the location of puncture (±). The parameters con¬ 
trol the steepness of the attenuation. We take the small¬ 
est possible power index p = 4 to achieve convergence of 
the solutions to the constraints. 

The second attenuation function is used to modify the 
background mean curvature and the source term in the 
momentum constraint equations. This takes the form 


K = /(+)si^(+) -b f(-)gK(_) , 


where 


9 

9± 


G{r±) 


5-t X 5- , 


{ 1 if r± > Tinax 

0 if r± < r^in , 

G{r±) otherwise, 


1 

2 


1 -b tanh 



TT 

2 


-1 


r± - r„ 





and the parameters < Tmax are chosen to be within 
the horizon. Note that the g attenuation function is used 
to modify the constraint equations themselves inside the 
horizon. We will refer to the above data as the standard 
data in the sections below. 


In addition, we consider a second type of initial data 
closely related to the above approach. For this, which 
we shall refer to as approximate data in the sections be¬ 
low, we analytically remove the singularity associated 
with DiM^^ by making the following two approximations. 
First, we take Mij to be the sum of the two Kerr trace- 
free extrinsic curvatures without correcting for the fact 
that the background metric is now a superimposed met¬ 
ric. Second, in the source term of Eq. (11b), we replace 
D,A*i with (A - Df)A';l + where Df and 

are the covariant derivative and extrinsic curvature asso¬ 
ciated with the two background conformal Kerr metrics. 
The former term contains no derivatives of A ’’^, while the 
latter is evaluated analytically. This is an additional ap¬ 
proximation because we neglect the fact that indices in 
are actually raised using the full superposed metric, 
rather than the associated Kerr metric. The net result 
of these approximations is that the initial data solution 
converges (relatively) rapidly with collocation points, but 
the resulting constraints converge to a small non-zero 
value. We use a subsequent CCZ4 evolution to remove 
this residual violation. This allow us to quantify the ef¬ 
fects of small violations of the initial constraints and how 
to control them. 


As shown in Fig. we verified the exponential con¬ 
vergence of the constraint violations of the initial data 
(with collocation points) using the L^-norms (RMS). We 
find that volume-averaged constraint violation (i.e., 
over the entire simulations domain) converge to levels of 
~ 5 X 10“^°. The constraint violations are largest near 
the x-axis (see Fig. [^. To verify the convergence of the 
data there, we calculated the L^-norm over as small box 
of width 0.5M centered on the x-axis. This box is cho¬ 
sen so that it lies just outside one of the horizons. The 
momentum constraints converge exponentially, but at a 
relatively slow rate, in this volume. The Hamiltonian 
converges exponentially to a level of ~ 10“^. The source 
of the relatively large violations is a high-frequency com¬ 
ponent in the initial data induced by the scale of the 
attenuation function. Note that the convergence of the 
approximate data is much faster with collocation points, 
but also converges to a non-zero value. 

Figure shows the initial data Hamiltonian con¬ 
straint violation on the xy-plane for a configuration with 
P/mirr = 1 and initial separation d/M = 10 . Note the 
high-frequency residual. By increasing the width of the 
attenuation function g above, we were able to partially 
mitigate the high-frequency noise in the constraint resid¬ 
uals using N = 192^ x 4 collocation points. For compari¬ 
son we also show the same configuration with P/mi„ = 1 
and initial separation d/M = 10 but for N = 48^ x 4 
collocation points, which represents a medium resolution 
for BY data. The choice of a lower number of collocation 
points {N = 48^ x 4) for BY is because the BY system 
is algebraically simpler than the new data (among oth¬ 
ers, the momentum constraints are solved exactly, and 
the background is flat). We thus expect that for a given 
number of collocation points, BY data will have a much 
smaller constraint violation, which is indeed what we see 
(see bottom panels of this figure). From the figure, we 
see that we can reach acceptable levels of constraint vio¬ 
lations with our data, but require a much larger number 
of collocation points than for BY. 


C. Evolution 

We evolve black hole binary initial data sets using the 
LazEv [28j implementation of the moving punctures ap¬ 
proach for both the BSSNOK formalism [29lj^ and the 
conformal and covariant formulation of the Z4 (CCZ4) 
system (Ref. [H]) which includes stronger damping of 
the constraint violations than the BSSNOK system. For 
the runs presented here, we use centered, eighth-order ac¬ 
curate finite differencing in space [32] and a fourth-order 
Runge-Kutta time integrator. Our code uses the Cac- 
tus/EinsteinToolkit |33l[34l infrastructure. We use 
the Carpet mesh refinement driver to provide a “moving 
boxes” style of mesh refinement [3^. Fifth order Kreiss- 
Oliger dissipation is added to evolved variables with dis- 
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FIG. 1. The L^-norms of the Hamiltonian and momentum constraints for the standard data over the full grid extending to 
400M (top-left), over a small volume centered on (x,y,z) = (5.0,0.35,0.75) (top-right), over a small volume containing the 
i-axis centered on (4.5, 0,0) (bottom-left). The number of collocation points is given by x x 32 in the A, B, and (j> 
dimensions, respectively. The bottom-right plot shows the constraint violation over the full grid with the approximate data 
instead. 


sipation coefficient e = 0.1. Note that when using CCZ4, 
we chose damping parameters ki = 0 . 1 , K 2 = 0 , and 
K 3 = 0 (see [H]). 

We locate the apparent horizons using the AHFind- 
erDirect code [3^] and measure the horizon spins us¬ 
ing the isolated horizon (IH) algorithm [37] . To compute 
the radiated angular momentum components, we use for¬ 
mulas based on “flux-linkages” |38| . explicitly written in 
terms of ^'4 jssmoj. We then extrapolate those extrac¬ 
tions to an infinite observer location using formulae ac¬ 
curate to Cl(l/rQj^g) [30] . 

We obtain accurate, convergent waveforms and horizon 
parameters by evolving this system in conjunction with a 
modified 1-f log lapse and a modified Gamma-driver shift 
condition [suniiii]. The lapse and shift are evolved with 

{dt - /3*a,)a = -a^f{a)K , (12a) 

- 7;/3“ . ( 12 b) 

where 77 = 2 . 

We have found that the choice /(a) = 8/(3q;(3 — 

a)) (approximate shock avoiding [ 21 ]) proves to be 

more stable and convenient when dealing with highly 
boosted moving punctures at relatively short separations, 
« lOOM (this proved particularly useful for the CCZ4 


simulations described below). This is due to the fact 
that the shock avoiding gauge suppresses a large am¬ 
plitude gauge wave that would otherwise be focused by 
the black holes and subsequently trigger a Courant vio¬ 
lation when the lapse gets too big. For the initial form 
of the lapse we use a{t = 0 ) = 1/(27/!bl — 1 ), where 
V’BL = 1 + m(+)/(2r(+)) -I- m(_)/(2r(_)). This proved 
to produce more accurate evolutions for highly spinning 
black holes m and we will also adopt it for the highly 
boosted cases in this paper. 

For both sets of evolutions, CCZ4 and BSSNOK, there 
are between 11 and 13 levels of mesh refinement depend¬ 
ing on the momentum of the black holes. Since we start 
the initial separations of the CCZ4 simulations much far¬ 
ther apart than the BSSNOK evolutions, the coarsest 
levels of the grid structure differ between the two sets of 
evolution. The CCZ4 evolutions have an outer boundary 
of 800M, while the BSSNOK evolutions have an outer 
boundary of 400M. The finest levels are the same for 
both sets of evolutions. We label the different resolution 
runs by nX where A is a global grid factor. For all evolu¬ 
tion runs, we use a grid of nl20. For this case, the finest 
resolution for the P/m\„ = 0.3 case is MjWl.2 and finest 
resolution for the Plm\„ = 4.0 case is Full 

details of the nl 20 grid structure is given in Table [T] 

Figure [^ shows the constraint violations versus time 
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FIG. 2. The Hamiltonian constraint in the a:y-plane of the standard data for Pjmirr = 1 with initial separation d/M = 10 in 
the region around the black holes for N = 192 collocation points (Top left). For comparison we also show the violations of the 
Hamiltinian constraint for BY data for N = 48 collocation points (Top right). In the lower panels we show the violation of 
the Hamiltonian constraint along the axis containing the black holes and perpendicular to it for our data and for BY data for 
different number of colocation points N = 32,48, 64. 


for a PjmiYr = 2 simulation using BSSNOK evolutions 
of stondarrf (fata for three resolutions (nlOO, nl20, nl44). 
During most of the run, after the initial settling of gauges, 
the constraint violations are convergent. We observe a 
hyperconvergent (4th-8th) order for the merger phase 
and then an slower convergence (due to residual grid in¬ 
terboundary radiation) of nearly hrst order after merger. 


III. RELATIVISTIC HEAD-ON COLLISIONS 

A major difference between our work here and pre¬ 
vious studies, see Refs. [101 dHHSllle], is that we use 
non-conformally flat initial data. The Bowen-York ini¬ 
tial data used previously are limited to representing black 
holes moving at speeds v < 0.9c, as shown in Fig.[^ The 
reason for this is that the assumption of conformal flat¬ 
ness introduces a Brill wave that gets stronger as the 
momentum parameter is increases. Most of this wave is 


absorbed by the black holes, leading them to increase in 
mass proportional to the momentum parameter. The net 
effect is that the ratio of momentum to mass of each black 
hole can never be larger than Plm\„ ^ 2 (i.e., vjc ^ 0.9). 
Note that here we use the irreducible mass of each black 
hole in place of the particle rest mass. 

The situation is similar to that observed in highly spin¬ 
ning black holes, where the conformally flat ansatz for the 
3-metric leads to a limitation [43H45j in the maximum in¬ 
trinsic spin of the black hole of around Sjw? « 0.93. 

On the other hand, Fig. shows that the new data we 
use here is not limited by this condition and can reach 
velocities closer to the speed of light, i.e. v ~ 0.99c. 
This is due to the much lower initial radiation content of 
the data. We will exploit this characteristic of the initial 
data to obtain a more accurate estimate of the maximum 
gravitational radiation produced by head-on collision of 
two equal mass, nonspinning, black holes. 

In order to explore the dependence of the radiated en- 












6 


TABLE I. Table of grid structure for case nl20. For P/rrn„ 
up to 2 we use up to mesh level 10. For P/rriirr = 3 we 
include an additional level and for P/mirr = 4 we use all 13 
mesh refinement levels. Whether the refinement level’s grid 
is centered on the origin or around the black holes (BHs) is 
given in column 2. The radius of the box is given in column 
3. For meshes with two values, the first is for the BSSNOK 
evolution of the standard data, and the second is for the CCZ4 
evolutions of the approximate data. 


Mesh Number 

Gentered on 

Radius 

Resolution 

0 

Origin 

400,800 

M/0.3 

1 

Origin 

200,500 

M/0.6 

2 

Origin 

140,300 

M/1.2 

3 

BHs 

32 

M/2.4 

4 

BHs 

16 

M/4.8 

5 

BHs 

8 

M/9.6 

6 

BHs 

4 

M/19.2 

7 

BHs 

2 

M/38.4 

8 

BHs 

1.2 

M/76.8 

9 

BHs 

0.6 

M/153.6 

10 

BHs 

0.3 

M/307.2 

11 

BHs 

0.15 

M/614.4 

12 

BHs 

0.08 

M/1228.8 



m 

FIG. 3. The L^-norm of the Hamiltonian constraint violation 
versus time for a P/mirr = 2 simulations using standard data 
evolved with BSSNOK at three different resolutions. At early 
times, the constraint violations are much smaller, but the 
violations become much larger than their initial values for 
most of the simulation. Globally the constraints decrease with 
resolution only slightly due to high-frequency noise on the 
grid. 


ergy on the magnitude of initial momentum of the two 
black holes, and then extrapolate the results to the ultra- 
relativistic limit, we set up a series of simulations with 
P/miir ranging from 0.3 to 4.03 (see Tables [n| and fill. 
We have chosen a relatively large initial separation of the 
black holes in order to ensure that the isolated horizon 
formalism can be used to accurately measure the mass 
of the black holes and to ensure that the momentum pa¬ 
rameter used in the simulations corresponds closely to 



Y 

FIG. 4. The center of mass speeds of the black holes for 
the Lorentz boosted data in this paper and Bowen-York 
data. The latter displays a limitation, reaching speeds of only 
V < 0.9c, while the former can reach near the ultrarelativistic 
regime. The thin line represent the special relativistic speed 
expression v = c-y/l — l/y^. 


the momentum of the black holes at infinite separation. 

In Fig. we show the {£ = 2,m = 2) mode of ipA for 
a typical simulations (here P/mi„ = 2) as seen by an 
observer at r = 130M. The spurious radiation is evident 
in the burst near t ^ 150M, well before the merger sig¬ 
nal. As can be seen from the figure, for the new data, 
the spurious radiation contains about 4% of the total en¬ 
ergy radiated for the standard data. On the other hand, 
the Bowen-York spurious radiation content is 24% of the 
total. 

To extrapolate the energy radiated to infinite observer 
location, we use 7 finite observers and extrapolate using 
a 1st order and 2nd order polynomial. 

For the standard data simulations (which were all 
evolved using BSSNOK), the extraction radii extended to 
J'obs = 130M. The error in the extrapolation is estimated 
by the difference between the two fits and is labeled “Inf 
Radius” error in Fig. 

For the simulations of the approximate data (which 
were all evolved using CCZ4), we needed to use larger 
initial separations than for the standard data in order to 
reduce the constraint violations on the initial slice. We 
therefore extracted the radiation at correspondingly large 
distances. For example, for the largest separation run 
d = 400M, the largest extraction radius as 275M (note 
the black holes were initially located at a: = ±200M). 

In all of the CCZ4 cases, the extrapolation formula of 
Ref. [20] to ©(l/r^j^g) gives a very robust set of values 
for the radiated energy. To provide a generous bound, 
we used those two radii as estimates of the infinite radius 
energy radiated. 

The other source of error we seek to keep under control 
is the initial, unphysical radiation content. We checked 
this spurious radiation for all of the Lorentz boosted 
runs. To do this, we compare the radiated energy of the 
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TABLE II. Table of initial parameters and energy radiated for the standard initial data evolved with BSSNOK. 


P/Madm 

Maum/M 

mirr/MADM 

P jvri'YP'p 

7 

d/M 

Erad /MaOM 

^Erad/A/ADM 

0.1437 

1.0008 

0.4804 

0.30 

1.0438 

100 

0.0011 

4.8e-5 

0.2238 

1.0028 

0.4488 

0.50 

1.1174 

100 

0.0031 

1.3e-6 

0.3547 

1.0093 

0.3555 

1.00 

1.4126 

100 

0.0182 

2.6e-4 

0.4510 

1.0177 

0.2583 

2.00 

2.2336 

100 

0.0585 

1.3e-3 

0.4792 

1.0268 

0.1597 

3.00 

3.1630 

100 

0.0858 

1.8e-3 

0.4886 

1.0250 

0.1220 

4.00 

4.1272 

150 

0.0957 

1.3e-4 


TABLE III. Table of initial parameters and energy radiated for the approximate data evolved with CCZ4. For each system, the 
initial ADM mass is normalized to 1. 


P/Madm 

Madm/AI 

mirr/AlADM 

P /^irr 

7 

d/MADM 

Prad /MaDM 

5Prad/MADM 

0.1439 

1.0000 

0.4807 

0.30 

1.0440 

100 

0.0011 

6.8e-7 

0.2245 

1.0000 

0.4498 

0.50 

1.1180 

100 

0.0030 

6.1e-6 

0.3558 

1.0000 

0.3559 

1.00 

1.4142 

200 

0.0183 

1.2e-4 

0.4530 

1.0000 

0.2263 

2.00 

2.2361 

200 

0.0592 

4.7e-4 

0.4800 

1.0000 

0.1594 

3.01 

3.1717 

300 

0.0859 

l.le-3 

0.4908 

1.0000 

0.1217 

4.03 

4.1231 

400 

0.0988 

9.7e-4 


full waveform with that obtained by removing the initial 
transient. The effect of the initial transient is to change 
the total radiated energy by ~ 1.6% (relative to the total 
radiated energy). The effect of this spurious radiation 
on the accuracy of the total radiated energy is shown in 
Fig.|6] under the label “Spurious”. 

It is worth noting here that the waveforms are ex¬ 
tracted by a multipole decomposition at the observer lo¬ 
cation. In practice a few of the lower modes are neces¬ 
sary for an accurate account of the total radiation. For 
instance, the £-mode contributions to the CCZ4 simu¬ 
lations for a Pjmirr = 3 run (with initial separation 
d — lOOM) at robs = 275M gives that £ = 2 contains 
90%, £ = 4 contains 8.3%, and £ = 6 contains 1.68% of 
the total energy radiated. Thus our results will include 
modes up to £ = 6. 

To test the accuracy and consistency of our simula¬ 
tions, we performed a convergence study of the radiated 
energy, the main physical quantity studied here, for six 
runs (all with initial P/min- = 0.5). We increases the 
resolution in stepsizes of 1.2 between each run. The re¬ 
sults of evaluation of the final mass of the black hole from 
the measurement of the gravitational radiation losses is 
shown in Fig.[^ While the differences with resolution are 
small, they are compatible with the expected 4th order 
convergence of the evolution system. In this figure, we 
fit the data to the form y = ao + aix^, where oq and oi 
are fitting constants and p is taken to be 2, 4, and 6. 

Also for runs with initial P/mi„ = 2, the agreement 
between the radiated energies, as measured from the 
waveforms (extrapolated to observer location to infinity 
via [20]) and those inferred from the initial ADM mass 
minus the remnant horizon mass, provides a consistency 
measure for the numerical simulation (here we increased 
the resolution by factors of 1.1). Assuming the differences 
scales like ah^, the 6—power of convergence for the three 
highest resolution h runs is found to be 4.17697 ± 1.139. 


In addition, we fit the same data to the form ah^ + c, 
where b was fixed to 1, 2, 3, 4, and 5, and for each choice 
of b, we fit to a and c. The results are summarized in 
Fig.lH which shows different orders of extrapolation to 
infinite resolution. The best results are near the expected 
4th order convergence. Here, we find consistency in the 
final mass to within 5 x 10“^M. 

Another interesting aspect to explore is how appropri¬ 
ate the standard moving puncture gauges ( 12 ) are for 
evolving highly-boosted black holes. We found that at 
relatively short initial distances the BSSNOK formalism 
generates a gauge wave focused by the two black holes 
that then induces a large change in the lapse. This can 
drive the lapse beyond 0 = 1 and trigger a Courant 
violation. This problem was resolved by starting the 
black holes at larger initial separations, allowing the large 
gauge waves to sufficiently dissipate before the collision. 
We also found it beneficial to use an initial lapse of the 
form ao = 1 /( 2 '(/ibl ~ 1 ) and the approximate shock 
avoiding gauge profile /(a) = 8/(3a(3 — a)) (which we 
used for all CCZ4 simulations). 

When fitting the radiated energy as a function of the 
initial momentum we assume a relative error of 1 % in 
each computed energy to determine its weight in the fit. 
We fit these energy values as a function of the variable 
where 771;^,. stands for the irreducible mass of each 
initially boosted black hole with momentum ±P. The 
upper panels in Fig. and IT display the results of the 
fitting, assuming the dependence of the energy radiated 
is given by the ZFL behavior [13 as] 


— -F A + V (l-47^)log(7+\/7^-l) \ 

M °° 272 + 273^^^ ) ’ 

(13) 

where 7 = -^ 1 - 1 - {P/m[„)‘^ and the only fitting parame¬ 
ter is Eaa- The relative deviations are mostly below 10%, 
and in particular are around 2 % for the most energetic 



























P/tTlir, 


FIG. 5. Top panel: The {£ = 2,m = 2) mode of ^4 for the 
Pjmirr = 2 Bowen-York and Lorentz-Boosted simulations. 
Note the spurious signal at t ~ 150M. The inset zooms in on 
the spurious signal of the Lorentz-Boosted waveform. Bottom 
Panel: The energy radiated for the two waveforms in the top 
panel. The contribution of the spurious radiation can be seen 
by looking at the energy radiated up to time t ~ 200 M. 


simulated collision. 

To assess the dependence with the chosen fitting func¬ 
tion, we have assumed a fit of the form (y = Aexp[—i? x]) 
with two fitting parameters {A and B), y and x being the 
independent and dependent variables, i.e. i^’radZ-^ADM 
and mirrfP, respectively. The results of this fit are 
displayed in the lower panels of Fig. and 10 In 
spite of introducing two fitting parameters, we observe 
that the residuals are larger than the fit using the ZFL 
form (13), thus rendering further support to this behav¬ 
ior. We have also experimented with fittings of the form 
{y = Aexp[—Bx^]), introducing a third parameter C in 
the fitting function, and also assuming (7 = 2, but none 
of these options displayed better behavior than the ZFL 
choice. 

In either case of the fits shown in Fig. 10 the esti¬ 
mated maximum radiated energy is around 13%, which 
provides a robust estimate, all errors considered, of the 
form F^max/ATADM = 0.13 ± 0.01. 


FIG. 6 . Estimated errors for each component of the compu¬ 
tation of the radiated energy. “Inf Radius” is the extrapola¬ 
tion from the finite extraction radius robs to infinity. “Spuri¬ 
ous” is the effect of the initial radiation content of the data. 
“Truncation” is an estimate of the finite difference resolution 
used in the simulation, and “AH vs. GW” is a consistency 
measure of the radiated energy as computed by the gravita¬ 
tional waveforms or the remnant mass of the final black hole. 
Shown in cyan is the total energy radiated for that simula¬ 
tion. Top panel is the BSSNOK evolutions of the standard 
data, and bottom panel is for the GCZ4 evolutions of the ap¬ 
proximate data. 


IV. CONCLUSIONS AND DISCUSSION 


Using improved full numerical techniques, we have 
been able to provide a more accurate determination of the 
maximum gravitational radiation produced in the head- 
on collision of nonspinning black holes. These techniques 
utilize initial data for highly boosted black holes m with 
much less radiation content than the Bowen-York coun¬ 
terparts, and reach near the ultrarelativistic regime with 
speeds much closer to c. We have successfully extrap¬ 
olated the extracted waveforms to infinite observer lo¬ 
cations with the techniques of Ref. [20], and added up 
to £ = 6 modes in the computation of the radiated en¬ 
ergy. The evolutions of the initial data have been carried 

























































































































9 



8x 

FIG. 7. The convergence for the Pjmiii = 0.50 final mass 
calculated from the gravitational modes {tmax = 6). sixth- 
, fourth-, and second-order fits to the data are shown. The 
fourth-order fit is closest to the data. For this convergence 
study, we used six runs (n72, n80, n90, nlOO, nl20, nl44). 


mass difference 


0.001 o; 
0.0008 t 
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0.00041 
0.0002 
0.0000 
- 0.0002 
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FIG. 8. The difference between the final horizon mass as 
calculated using the IH formalism and as inferred from the 
radiated energy for a P/mirr = 2 simulation as a function 
of resolution. The three highest resolution runs are shown. 
The data points themselves appear to be convergent. We es¬ 
timate the infinite resolution limit by assuming second, third, 
fourth-order, and fifth-order convergence. The agreement be¬ 
tween the horizon derived mass and radiation inferred mass 
at infinite resolution is 5 x 10“®M for the expected 4th order 
convergence. 


out using the moving punctures approach using both the 
BSSNOK and CCZ4 systems. 

We find a maximum radiated energy of 13 ± 1% of the 
total mass of the system, with most of the errors coming 
from the functional fitting and subsequent extrapolation 
to infinite boost. This result is in close agreement 
with the analytic estimates of 13.4% of Ref. [T7] using 
thermodynamic arguments and the previous numerical 
estimate of 14 ± 3% in Ref. m- However, they seem 
to be in conflict with the analytic estimates of 16.4% 
from second order perturbations [1] and 17% from the 


multipolar analysis of the ZFL m- 




tTljrr/P 


FIG. 9. Fits to the energy radiated at infinity for the BSS¬ 
NOK evolutions of the standard data in Table [II] Upper 
plot: Fit using the 1-parameter ZFL like fit. Lower plot: 
An alternative two-parameter {A and B) fit of the form 
{y = Aexp[—Pa:]). Both fits use data assuming a weighting 
error of the points of 1% and include fits to both the energy 
radiated as measured by extraction of radiation (WF) or by 
the remnant mass (AH). 
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FIG. 10. Fits to the energy radiated at infinity for the 
CCZ4 evolutions of the approximate data in Table Iml Up¬ 
per plot: Fit using the 1-parameter ZFL like fit. Lower 
plot: An alternative two-parameter {A and B) fit of the form 
{y = Aexp[—Ba;]). Both fits use data assuming a weighting 
error of the points of 1% and include fits to both the energy 
radiated as measured by extraction of radiation (WF) or by 
the remnant mass (AH). 
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